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Abstract 



Left-right symmetry breaking is critical to vertebrate embryonic development; in many species this process begins 
with cilia-driven flow in a structure termed the 'node'. Primary 'whirling' cilia, tilted towards the posterior, trans- 
port morphogen-containing vesicles towards the left, initiating left-right asymmetric development. We review recent 
theoretical models based on the point-force stokeslet and point-torque rotlet singularities, explaining how rotation and 

Q^ surface-tilt produce directional flow. Analysis of image singularity systems enforcing the no-slip condition shows how 

tilted rotation produces a far- field 'stresslet' directional flow, and how time-dependent point-force and time-independent 
point-torque models are in this respect equivalent. Associated slender body theory analysis is reviewed; this approach 
enables efficient and accurate simulation of three-dimensional time-dependent flow, time-dependence being essential in 
predicting features of the flow such as chaotic advection, which have subsequently been determined experimentally. A 
new model for the nodal flow utilising the regularized stokeslet method is developed, to model the effect of the overlying 
Reichert's membrane. Velocity fields and particle paths within the enclosed domain are computed and compared with 
the flow profiles predicted by previous 'membrane-less' models. Computations confirm that the presence of the mem- 
brane produces flow-reversal in the upper region, but no continuous region of reverse flow close to the epithelium. The 
stresslet far-field is no longer evident in the membrane model, due to the depth of the cavity being of similar magnitude 
to the cilium length. Simulations predict that vesicles released within one cilium length of the epithelium are generally 
transported to the left via a 'loopy drift' motion, sometimes involving highly unpredictable detours around leftward 

• • cilia. Particles released just above the cilia tips were not predicted to reach to the extreme edges of the node, but rather 

^ are returned to the right by the counterflow. Flow to the right and left of the cilia array is of very small magnitude, 

KJi suggesting that effective transport of particles to the extremities of the node requires cilia to be distributed all the 

^ Y^ way to the edges. There is no continuous layer of rightward flow close to the epithelium, except for a region close to 

^ the posterior edge of the node. Future work will involve investigating issues such as the precise shape of the node and 

cilia distribution and the effect of advection and diffusion on morphogens, hence explaining more fully the role of fluid 
mechanics in this vital developmental process. 



1 Personal reflection 

It is a great honour to be invited to contribute to the volume in memory of Professor Ernie Tuck, the last authors' 
undergraduate project supervisor, mentor and friend for life, his influence being evident in this paper. During the period 
1968-69 Tuck introduced the last author to slender body theory in Stokes flow, based on his earlier papers [1, 2 . Part of 
his project work has recently been published, 40 years on [3 . The application of slender body theory to micro-organism 
locomotion led to studying ciliary propulsion, the study of which began during Tuck's period at Caltech under the direction 
of the distinguished hydrodynamicist, T.Y. Wu. While at Caltech, Tuck [4 extended the famous Taylor [5 analysis of a 
swimming sheet in a viscous fluid by adding inertia and tangential motion to the analysis. As a result Tuck introduced the 
last author to the paper by Lighthill [6 on the squirming motion of a sphere in a viscous liquid, which was found to contain 
an analytical error. This later led to the last author undertaking a Ph.D under the supervision of Sir James Lighthill, 
the recently appointed Lucasian Professor of Mathematics at the University of Cambridge. One of the first actions, 
encouraged by Lighthill, was to publish a corrected version of the paper [7 but with an application directed at ciliary 
propulsion. Lighthill [8 reported on the outcome of some of these studies as one of the invited speakers at the ICTAM 
Congress in Moscow in 1972. Thirty six years later, in 2008, Ernie Tuck was the president of the organising committee of 
the highly successful ICTAM Congress held in Adelaide. This paper will incorporate the ideas originally initiated during 
this 1968-69 period in the study of one of the most exciting current areas of science, embryonic development, an area 
which has previously attracted the attention of such luminaries as Alan Turing [9 and Francis Crick [10] . 

2 Introduction 

We begin by reviewing the relevant biology of nodal cilia and left-right symmetry-breaking, and briefly outlining experi- 
mental and theoretical research into this process. In subsequent sections of the paper we review current knowledge of the 
fluid mechanics in more detail, focusing on the use of singularities in explaining and modelling the flow. We then develop a 
new model which takes into account the effect of an overlying membrane, and present computational calculations showing 
its effect on the fluid flow and particle transport. 

On first impressions, the bodies of mammals are essentially bilaterally symmetric. Beneath the surface it is very 
different however, with symmetry being broken in an organised way, the human heart for example being on the individual's 
left, the liver on the right (figure la). The conventional left-right axis designation in this figure, with the 'left' of the 
individual appearing on the right of the figure will be used throughout this paper. This symmetry and asymmetry is laid 
down in the early stages of development, and is critical to embryonic survival. In this section we shall briefly describe the 
relevant biology, before addressing the fluid mechanics; for detailed review, see Hirokawa et al. [11 . 

The breaking of the left-right axis occurs after the establishment of the anterior-posterior and dorsal-ventral axes, 
but until the mid 1990s there was little understanding of the mechanisms involved. Kartagener in the 1930s, and others 
previously (for historical review, see [15, 16 ), had identified a triad syndrome of respiratory problems, male infertility and 
situs inversus, the latter being the lateral transposition of the internal organs (figure lb). The common factor causing 
these symptoms was not however understood. A serendipitous discovery occurred in 1974 when Afzelius examined the 
doctoral thesis of H. Pedersen. A link was made between Kartagener's syndrome and defects in the motor protein dynein. 
This defect results in immotility of cilia, such as those in the respiratory system, and of sperm flagella, which are required 
for motility and fertilisation. 

However, the connection between cilia and asymmetric development was not fully established for another twenty years. 
In 1994, Sulik et al. [17 reported on the development of a 'node' structure on the ventral side of mouse embryos at 7-9 
days post fertilisation (figure Ic, marked 'VN'), on which primary cilia were expressed (figure Id, e). Primary cilia are 
microscopic hair-like organelles found on virtually all cells in the body; they typically occur one-per-cell, their internal 
structure differing from that of mucus-propelling lung cilia, and sperm flagella, in that the central microtubule pair are 
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Figure 1: The role of cilia in generating the mammalian body plan, (a) Normal organ placement and (b) mirrored organ 
placement in situs inversus (copyright figure not included in arxiv version), (c) The ventral side of the mouse embryo at 
eight days post-fertilisation, with Reichert's membrane removed, showing the anterior-posterior and left-right axes, the 
ventral node being the depression structure highlighted (copyright figure not included in arxiv version). Abbreviations: 
VN, ventral node; NP, notochordal plate; FG, foregut. (d) Schematic of the ventral node, showing how the biological axis 
designations correspond to the Cartesian axes used in the mathematical model of this paper. The node is shown looking 
down onto the ventral surface of the embryo; the cones represent the approximate beat pattern of the cilia, (e) Electron 
micrograph of the nodal fioor and nodal cilia (copyright figure not included in arxiv version), (f) Schematic section of the 
node cavity, viewed from the posterior, looking towards the anterior, with the 'left' of the embryo being conventionally 
shown on the right of the diagram. The 'leftward' mainstream fiow is balanced by a rightward return fiow induced by 
Reichert's membrane; the existence of a rightward return fiow close to the node fioor is a subject of continuing discussion 



absent (detailed in [H]). 

In contrast to mucus-propelling cilia, primary cilia were generally thought to be non-motile. Sulik et al. hypothesised 
that motility of the primary cilia in the node may be required for left-right symmetry-breaking, this being the reason 
that situs inversus occurs in patients with defective dynein. Nonaka et al. |18] subsequently revealed that these cilia 
were indeed motile, and in normal embryos perform a 'whirling' motion that creates a right-to- left flow (figure Id, f). 
Gene knockout mouse embryos were produced which did not assemble nodal cilia; these mutants did not break left-right 
symmetry correctly and did not survive. 

Many significant advances have been reported in the biological literature since; we restrict our attention to four 
experimental and four theoretical studies relevant to the fluid mechanics. Nonaka et al. [19] showed that artificial flow 
could reverse the situs of normal embryos, and direct the situs of mutants with immotile cilia. A number of questions were 
posed: firstly, how do whirling cilia create such a directional flow, and secondly, how is the flow converted to asymmetric 
development? 

An answer to the first question was proposed by Cartwright et al. [20 . A fluid dynamic model using 'rotlet' point- 
torque singular solutions of the Stokes flow equations was formulated, from which it was deduced that a posterior 'tilt' 
in the rotational axis and the combined effect of multiple cilia arranged in an array could produce a right-to-left flow. 
Although this model did not take into account all of the essential physics of the cilium/surface interaction [21 , the 
'posterior tilt' was nevertheless a fundamental theoretical prediction, subsequently confirmed experimentally [22 (for 
more discussion of this work and the underlying modelling philosophy, see [23]). Nonaka et al. [13] devised a different 
approach to investigate the fluid mechanics of nodal cilia: a laboratory analogue model using tilted rotating wires, with 
fluid viscosity sufficiently high and rotational speed sufficiently low to produce a flow with appropriately small Reynolds 
number (figure |2k). A rotating rod tracing out a conical shape, tilted towards the posterior, can be defined with two 
parameters, the tilt angle, denoted ^, and 'semi-cone' angle, denoted ip (figure^). Nonaka et al. observed that tilted 
rotation in their wire cilia model did indeed result in directional transport of marker particles; moreover they investigated 
how the tilt and semi-cone angles affected the particle transport rate. Optimal transport occurred when the angle sum 
was close to 90°, corresponding to the cilium being close to the epithelium during the return stroke; moreover a semi-cone 
angle of 60° resulted in more rapid transport than semi-cone angle of 45°. The question of optimality of the cilia motion 
is discussed further in section [44l 

The second question, regarding the conversion of the nodal flow to asymmetric development, was subject to two 
different hypotheses, one involving morphogen protein transport, the other mechanical sensing by cilia. An important 
breakthrough was made in the same year by Tanaka et al. [24 , who found that lipid-enclosed parcels of morphogen 
proteins, that affect cell differentiation and development, are released from the nodal cells and pushed into the flow field 
by microvilli. These 'nodal vesicular parcels' (NVPs, figure If) were hypothesised to act as vessels for the transport of 
morphogens selectively to the left side of the embryo. The transport and break-up of these parcels within an enclosed 
volume was subsequently modelled as a steady flow in an enclosed domain, and solved numerically by the finite element 
method by Cartwright et al. [25j . 

Our first study [26 was designed to address two important physical effects: the role of cilium-cell surface interaction, 
and the effect of unsteadiness in the flow field. The methodology, reviewed in this paper, was based on a numerical 
application of slender body theory to calculate time- varying solutions to the Stokes flow equations, modelling the viscous- 
dominated microscale domain. While these equations have no explicit time-dependence, the whirling motion of the cilia 
does have an explicit time-dependence, producing an unsteady flow. An alternative explanation for nodal flow generation 
has been given in the literature, involving velocity differences in the effective and recovery stroke. However, due to the 
lack of explicit time-dependence of the Stokes flow equations, a motion performed in opposite directions with different 
speeds will not produce any net transport [27 , a fact sometimes referred to as the 'Scallop Theorem'. By contrast, the 
unsteady boundary conditions and resulting unsteady flow have important implications for the transport of particles 
within the node. 

Our computational simulations [26 investigated the flow produced by a small array of tilted cilia projecting from 





Figure 2: (a) Schematic figure showing the experimental wire ciha model, redrawn from Nonaka et al. [13]. (b) Schematic 
showing the coordinate system used for the mathematical model of the conical rotation, its relation to the biological axes, 
and the conical rotation model for an individual cilium. The angle is the posterior tilt, ?/^ is the semi-cone angle. 



a plane in a semi-infinite fiuid. This allowed calculation of three-dimensional particle transport under the unsteady 
fiow field, and showed that 'chaotic' effects could occur, specifically that particles close to the cilia beat envelopes, that 
are initially close together, could diverge considerably. This prediction was subsequently confirmed experimentally in 
zebrafish embryos using a novel laser ablation technique ( [28] 



4.5 for more detail). The results showed that 



see section 

an ensemble effect is not required to produce a directional fiow; a single cilium can produce a directional fiow through an 
image-induced stresslet far-field, as described in detail in section [42] The approach of using singular solutions and slender 
body theory was developed further in [21], wherein an analytical formula was developed for the volume fiow rate, closely 



paralleling the experimental findings of Nonaka et al. p^j this is described in section 4.4 Unsteady simulations of fiow 
due to relatively large arrays of cilia were also presented; however a limitation of this work is that it did not account for 
the effect of the upper Reichert's membrane (figure If) on the fiow field. This had previously been accounted for in the 
steady fiow model of Cartwright et al. [25], who reported two separate regions of return fiow: an upper return fiow some 
distance above the cilia tips (figure If), and an additional lower region close to the epithelium (not shown), driven by the 
return stroke of the cilia. It therefore remains to investigate the effect of the membrane in the context of a time-dependent 
model. We shall first review in detail the anatomy and geometry of the node, then describe the mathematical modelling, 
and finally present new results for time-dependent fiow with the upper membrane taken into account. 



3 Anatomy of the node and geometric model 

The most commonly-studied species is the mouse; the essential details apply to organising structures found in the majority 
of species studied. Again, for detailed review, we refer to Hirokawa et al. [11]. Briefiy, the node of the mouse is 
an approximately triangular depression that appears on the ventral side with the triangle apex pointing towards the 
anterior (figure lb, c), 50-100 fiva in width and 10-20 fiva in depth, covered with Reichert's membrane, and filled with 
extraembryonic fiuid. The ventral surface, forming a tissue termed 'epithelium', has a few hundred 'nodal pit cells', 
generally expressing a single cilium of length 3-5 fivn and diameter 0.3 /im, giving a slenderness ratio 0(1/10). The 
characteristic rotation rate has been quoted as approximately 600 rpm [11^, or 10 Hz. The 'nodal fiow' first reported by 



Figure 3: Experimental data of cilia distribution and kinematics in three species, mouse, rabbit and medakafish (copyright 
figure not included in arxiv version) (a) Schematic of the tilt angle ^, semi-cone angle ^Z^, angle of cilium rotation axis 
from anterior ^, cilia projected length p and rotation angle ^ — the latter three scalar variables are used only for this 
figure and are not explicitly referred to in the mathematical model. Axis directions are marked A, P (anterior-posterior) 
and L, R (left-right), (b) Data for tilt angle 6 versus angle of rotation axis from posterior, (j) — the approximate mean 
of 180° corresponds to a complete posterior tilt, (c) Frequency distribution of cilium projected length, (d) Frequency 
distribution of semi-cone angle, inferred from the shape traced out by the cilium tip during the beat cycle, (e) Data 
showing the rotation angle traced out over time by a population of cilia in the mouse node, with the mean rotation over 
one beat cycle shown. 

Nonaka et al. [18 was visualised as the leftward movement of tracer particles at 2-5 /im/s at around 5 jiia above the cell 
surface; a slower 'return flow' occurs closer to the overlying membrane ([20l[22]; figure If). 

To model the flow mathematically, we represent the epithelium surface by the plane X3 = 0, the cilia extending into 
the region X3 > 0, with the positive Xi direction corresponding to the 'left' and the positive X2 direction corresponding 
to the anterior. The cilia are therefore tilted towards the negative X2 direction. The cilia motion is modelled as a rod 
moving through a conical envelope; the cone has semi-cone angle ip and tilt angle 6. The cilium cannot penetrate the 
epithelium surface, hence ?/^ + ^ < 90°. This is summarised in figure |2]3. Experimental data reported by Okada et al. 
[22] of the variation in the angles and ^Z^, and additionally angle of the tilt axis from the anterior, denoted (/), are given 
in figure [sja, b, d) for mouse, rabbit and medakafish. Additionally, values of the apparent cilium length p are given 
(figure Isfc), showing considerable variation even within the class mammalia. The rotation rate, typically approximately 
10 Hz (figure |3^), varies by approximately 20% over a beat cycle in the mouse, slowing slightly during the recovery stroke, 
during which surface-induced viscous drag will be greatest. 

In the sections below we describe the mathematical theory, starting with the fundamental singularities of Stokes flow, 
and the image singularities associated with enforcing the no-slip condition on the cell surface. We then describe the use 
of slender body theory, and review the derivation of an analytic formula estimating volume flow rate due to a single 
cilium. We then discuss the physical aspects arising from the image systems and slender body theory in more detail, and 
how their predictions have subsequently been confirmed by experimental observations. We then apply the regularized 
stokeslet method to model Reichert's membrane for the first time, and report computational results showing how this 
modifies the flow fields and particle tracks. We finish by relating these simulation results to the biology of left-right 
symmetry-breaking, and outline future work for fluid mechanics and morphogen diffusion studies. 

4 Singularity methods for fluid mechanics modeUing 

The mathematical techniques in this paper will be based on singular solutions of the zero Reynolds number fluid flow 
equations, and their regularized counterparts. Firstly we review these solutions and how they give insight into the nodal 
flow, before describing the associated slender body theory. 

4.1 Fundamental singularities of Stokes flow 

Due to the very small length and velocity scales, the Reynolds number of the nodal flow is much smaller than one, and 
so the Stokes flow equations are a very accurate model for the fluid dynamics. In this section we review a number of 
'fundamental solutions', representing the flow due to concentrated force, torque, strain and a source dipole. The linearity 
of the Stokes flow equations allow these solutions to be summed or integrated to satisfy boundary conditions associated 



with the epithehum or Reichert's membrane, and the flow produced by the movement of the ciha. Hence the properties 
of these solutions may be used to understand the physical mechanisms generating the flow, and also as a basis for slender 
body theory and boundary integral methods, which allow efficient and accurate computational modelling. 
The Stokes flow equations with a force per unit volume F centered at y are, 

Vj9 = /iV^n + Fx(cc-?/), V-n = 0, (1) 

where p = p{x^t) is pressure, u = u{x^t) is velocity and ji is dynamic viscosity. The term x{'^ ~ v) ^^Y be the Dirac 
delta distribution, representing a singular point-force, or a smoothed 'blob' function (section [5?T] ). With x(^ — u) taken 
as the delta distribution 5{x — i/), the solution to these equations is given by the stokeslet tensor, which in an inflnite 
fluid takes the form 

*'<-^> = 8i(T + 7?)- '^> 

Ti = Xi — i/i {i = 1, 2, 3) and r^ = rl-\-r2-\-rl; the velocity then being given by Ui = X].-^i SijFj. Henceforth the summation 
convention over repeated indices will be assumed. Taking x to be a smoothed blob function leads to a regularized stokeslet 
solution, introduced by [29 . The inflnite domain stokeslet decays with 0(l/r) in the far-fleld. 

As described in detail in Blake & Chwang [30^ , the stokeslet singularity may be differentiated to give a stokes-doublet 
singularity, a rank 3 tensor. This can be decomposed into symmetric and antisymmetric terms, referred to as the 'stresslet' 
and 'rotlet' singularities respectively 
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(3) 



The ffrst term is a rank 3 tensor, symmetric in j, /c, termed the 'stresslet'; when multiplied by a rank 2 tensor giving 
the stresslet strength, it generates a velocity ffeld due to a concentrated straining motion of the ffuid. The second 
term, antisymmetric in j, /c, termed the 'rotlet', can be shown to be a rank 2 tensor; when multipled by a vectorial 
strength, it generates the velocity field due to a concentrated torque applied to the fiuid. The rotlet can be rewritten as, 
Lij{x — y) = eijkrk/iSTTjur^)^ the velocity field being given by Ui = LijQj for point-torque Qj. Both the stresslet and 
rotlet velocity fields decay with O (l/r^) in the far-field, more rapidly than the stokeslet. 

Another relevant singular solution of the Stokes fiow equations is the point source dipole, a special case of a stokes- 
quadrupole, which tensorially is given by. 



A,(a.,y) = -^(^-^). (4) 



47r 

This can be derived by differentiating the stresslet, then multiplying by the idemfactor Sjk- The source dipole decays 
more rapidly than the stresslet or rotlet, being O (l/r^). This singular solution is used in both the construction of no-slip 



image systems in section 4.2, and the formulation of slender body theory (section [4. 3 [). 



4.2 From fundamental singularities to nodal flow 

An array of infinite domain rotlets Lij was used by Cartwright et al. ^20j, inspiring their prediction of the posterior tilt; it 
was shown that an array of rotlets produces a unidirectional fiow in the region spanned by the singularities. This model 
does not however satisfy the no-slip condition on the epithelium; moreover it predicts that the far-field fiow is rotational. 
In this section we review the work of Blake ^31j and Blake & Chwang [30 on the use of image singularities to satisfy the 
no-slip conditions, and the implications of this for the nodal fiow. 



To model cilia projecting from a cell surface, Blake [ST derived the following solution for a point-force a height h 
above the no-slip surface X3 = 0, 
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(5) 



where the tensor Ajk takes value +1 for j = /c = 1,2, value —1 for j = /c = 3, and zero if j ^ /c, originally written as 
{SjaSka — ^js^ks) with a = 1, 2 by Blake [31]. The image location is given by Ri = ri^ R2 = r2 and Rs = —h. 

The image system can be analysed as shown in figure [4J briefiy it consists of an equal and opposite image stokeslet, a 
stresslet and a source doublet. The far-field of Bij is that of a stresslet, and implies that the direction of the streamlines 
will asymptote to being parallel to the vector x — y. The streamlines of the stresslet point out radially in the far-field, 
which was evident in the particle tracking simulations of Smith et al. [21 . Moreover, the stokeslet decay rate of 0(l/r) 
is converted to O (l/r^) by the image system. The stokeslet and image system have been used for an array of problems 
both within and outside of biological fiuid mechanics; for the problem considered they form the basis of a slender body 
theory model of cilia protruding from a cell surface, as described in section [43) 

The rotlet singularity also has an associated image system, first given by Blake & Chwang [30] , 
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(6) 



Of note is that the far-field is similar to that of the stokeslet singularity, being a stresslet; a rotational near-field motion 
is 'converted' to a far-field straining motion by the infiuence of the surface. 

Stokeslets combined with their image system have generally been used for time-dependent models of cilia-driven nodal 



fiow using slender body theory, as in [26 . This approach will be described in more detail in section 4.3 The rotlet 
with no image system has been used as a basis for a time-averaged model [20]; from the above physical arguments we 
observe that once the rotlet image system is taken into account, the two approaches are consistent in the far-field fiow 
they produce. For further discussion of the consistency of these approaches, see Smith et al. [26 . 

4.3 Slender body theory 

Slender body theory represents the fiow driven by one or many cilia as centreline distributions of stokeslet and higher 
order singularities. For review and computational implementation, see [26j, for asymptotic justification, see [32 and for 
an analytic approach to the modelling of bodies with varying diameter axisymmetric geometries, see Blake et al. [3 . 
The essence of these methods is that the high aspect ratio of the cilium radius to length a/L allows both near-field and 
far-field fiows to be accurately approximated by a line distribution. 



Jo 



Ui{x,t)= / [Bij{x;^{s,t))fj{s,t)^Dij{x;^{s,t))gj{s,t)] ds, (7) 

Jo 

where fj{s^t) is force per unit length on the cilium, and gj{s,t) is a source-dipole distribution strength. The task is then 
to determine strengths of these line distributions that will satisify as closely as possible the boundary condition Ui{x) = 
dt(,i{s, t), where x is any point on the cilium surface around the centreline point ^(sq, t). The curved surface of the cilium 
can, for example, be represented conveniently as a curved cylinder, X(5o, t; a) := ^(sq, t) + a(n(so, t) cosa-\-b{so, t) sin a) 
where ^ a < 27r and n, h are unit normal and binormal respectively. In this case, the appropriate dipole weighting, 
as originally used by Hancock [33 is gj{s,t) = — (a^/4/i)/^(5, t), where fMs^t) is the component of the force density 
projected onto the plane given by the normal and binormal. This model is reasonably accurate away from the cilium tip 
(see [26l[27l Chap 6]); from it we can develop a simplified 'resistive force theory' model, as first used by Gray & Hancock 
[34] . A simple derivation follows: 
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Figure 4: Diagrams illustrating the image system for stokeslets (a) j = 1 and (6) j = 3, the strengths of the components 
being given in brackets. Note the significant difference in the far-fields; a stresslet for j = 1 or 2 and a lower order and 
weaker stokes-quadrupole for j = 3. (c, d) The image system and far-fields for a rotlet with components (c) tangential 
and {d) normal to the wall. 



Consider the velocity field on the cilium surface, and approximate the integral with its 'local' part only, 

-,2 



u 



:(X(5o,t;a)) = / 



\s-so\<q 



Bij{X{so,t;a);^{s))fj{s,t) + ^Dij{X{so,t;a);^{s))ff{s,t) 



ds, 



(8) 



for an appropriately-chosen length q. Neglecting cilium curvature on the segment {sq — <7, sq + ^), and assuming that the 
force density can be approximated by its midpoint value fj{so^t)^ the integrals can be calculated analytically, as first 
given by Hancock [33J, 



u,{X{so,t;a)) = (Ct)-' fUso,t) + {C„)-' fnso,t) + (Cb)-' fUso,t), 



(9) 



where the 'resistance coefficients' Ct^Cn^Cb in the tangential, normal and binormal directions respectively can be cal- 
culated analytically. It is not possible to derive optimal values of these in general [35 , however for cilium motion, 
successfully-used values are 

Q= . ..^:^. ,. , c„=c.= . ..f"i ,. , (10) 



-l + 21og(2g/a)^ 



l + 21og(2g/a)^ 



with q = {aLY''^ [36]. For flagellar motions and '9+2' cilia motion, an important property of these coeflicients is the 
anisotropy ratio 7 = Cn/Ct ~ 2. For our model of a nodal cilium as a straight whirling rod, tangential motion is zero, 
and hence the values of Ct and 7 are not relevant. This 'resistive force theory' approach will be used in section [44] to 
estimate the volume flow rate produced by a single cilium. The formulae for Cn and C5 show that the force density on 
the cilium is 0{ujL/\og{L/a)) where uo is the angular frequency of rotation. 

While this approach gives a good approximation to the cilium-fluid interaction, if accurate flow velocity flelds are 
required, this method is limited in that is not possible to obtain a force distribution that satisfles Ui{$,{sQ^t)) = dtii{sQ^t) 
uniformly towards the cilium ends. A more reflned approach is to use the quadratic dipole weighting gj{s,t) = — a^s(l — 
s)/j (s, t), based on the analytic solution of Chwang & Wu [37] for a straight rod, subsequently generalised to an asymptotic 
analysis by Johnson [32] and numerical implementation by Smith et al. [26]. If the surface of the cilium is taken to be a 
curved slender ellipsoid, with X(s, t; a) := ^{s^t) -\-a[l — {s — l/2)^/(a^ + l/4)]-'^/^(n(5, t) cos a + 6(5, t) sin a), the singular 

solutions being distributed between the foci of the ellipsoid, 5 = 0, 1, then accuracy to O l{a/L) ) can be obtained 
uniformly along the cilium, including the ends. In section [6] we denote the combination of plane boundary stokeslets Bij 
and source dipoles Dij with weighting gj{s^t) = — a^s(l — s)fj{s^t) by the symbol Gij. 

4.4 Volume flow rate 

To give an indication of the semi-cone and tilt angles that give optimal fluid transport, we consider the volume flow rate 
in the Xi ('leftward') direction across the half-plane given by X3 > and —00 < X2 < 00, for flxed arbitrary xi produced 
by a line distribution of the form (l7|). The source dipole Dij produces zero flow across this plane, while a unit stokeslet 
Bij located at X3 = /i pointing in the Xj direction produces a volume flow rate of /i/(7r/i) for j = 1 and zero otherwise 

EH]. 

The conical rotation can be parameterised as follows — with no posterior tilt. 



Ci(s,t) = 


- ssin'0cos(cL;t), 


Us,t) -- 


= — 5sin?/^sin(a;t)), 


Us,t) -- 


= 5 cos?/;. 



(11) 

Applying an anticlockwise rotation by angle d about the x\ axis, we have 

^i(s,t) = 5sin'0cos(cjt), 

^2{s^t) = s{—simljsm{ujt)cosO — cosiljsmO)^ 

^3(5, t) = s{—simljsm{ujt)smO-\-cosipcosO). (12) 

The velocity 1/(5, t) of the cilium can be determined by taking the time derivative of ^(5, t). Noting that 1/(5, t) is normal 
to the cilium centreline, we can deflne a unit normal vector n{t) such that u{s,t) = Un{s^t) = uosnit). Resistive force 
theory gives the approximation fi{s^t) = CnOJsni{t)^ with ni{t) = — simlj sm{ujt) . 

As discussed above, a point-force of strength /i pointing in the xi direction, situated at height ^3 above a no-slip 
surface, produces an instantaneous volume flow rate of /i^3/(7r/i); point-forces acting in the X2 and xs directions produce 
zero flow. The volume flow rate induced by the force CnOOsn{t) is therefore equal to Cni^sni{t)(,s{s,t)/{7rfi). Integrating 
this function from 5 = to L, and averaging over a beat cycle, the mean volume flow rate can therefore be written, 

Q = i^j J sm{t)Us,t)dtds. (13) 
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Using the above expressions for ^3(5, t) and ni{t) and noting that the integrals of sin^(co't) and sin(co't) cos(co't) over one 
period are 27r /uj and zero respectively, we find that the above integral evaluates to 

Q = -^ sin^ tjj sin 0, (14) 

67r/i 

as given in [2Tj. 

The mean volume fiow rate (3('0, 0) is the product of two non-negative, non-decreasing functions, and therefore attains 
a maximum value when -\- t/j = 90°, in agreement with the outcomes observed by Nonaka et al. [13 . Further analysis 
shows that the maximum occurs at '^ = arctanV^ = 54.7° and = 35.3°, ignoring the fact that the slender body model 
loses accuracy as the cilium closely approaches the surface. This occurs when -\- iIj = 90°, the surface drag in the 
real system becoming much larger than CnU^- A range of values of the tilt angle 6 have been observed experimentally 
[22] . however the mean values in mouse, rabbit and medakafish are typically between 35-40° (figure Isb). There is more 
spread in the reported values of the semi-cone angle ?/^ for these species, their means ranging from approximately 40-50° 
(figure §1). 

As discussed in the introduction, the experimental model of Nonaka et al. ^13^ using mechanically-driven wire 'cilia' 
(figure |2^) produced the most effective transport when 6> + ?/^ = 90°. Furthermore, comparing transport with semi-cone 
angles of ?/^ = 45° and ip = 60°, more effective transport occurred in the latter case, a result with which the functional 
form for Q is consistent [21 . 

We now review in more detail how slender body theory and fundamental solutions explain the physical mechanisms 
underlying fiow generation. 

4.5 Physics underlying the nodal cilium beat cycle 

The physical effects governing the fiuid motion generated by the whirling nodal cilia can be interpreted through equa- 
tion ([5| and figure |4ja, b), which give the velocity fields due to a force (stokeslet) near a rigid boundary. In figure [sja, 
b) the 'influence zones' are schematically illustrated for the nodal cilium in its effective stroke (upright orientation) and 
recovery stroke (near epithelial surface) based on the stokeslet model. The diagrams in figures [sja, b) are adaptations 
of earlier diagrams for the motion of planar-beating '9+2' cilia presented in Blake [39l|40], Blake and Sleigh [41] and for 
mucous transport by Sleigh et al. [42 . 

Figure [5^ shows the key physical features of the effective stroke highlighting the three zones of infiuence: 

• the inner-field — in the region for which r <C 1, the scaled cilium length, the velocity field is O(logr), as shown by 
integrating stokeslet and source-dipole distributions [43] , 

• the near-field — in the region for which r = 0(1), the velocity field is 0(l/r), essentially dominated by the stokeslet 
field, although as r becomes comparable with the distance from the image system 2/i, this will begin to have a 
significant effect. 

• the far-field — in the region for which r ^ 1, the stokeslet and image system are approximated by a symmetric 
stokes-dipole, and outcome of the equal and opposite force of the wall and the couple term in the image system, 
giving a 'stresslet' far-field, which decays with 0(l/r^). 

Figure JSJd shows the features of the recovery stroke, which is close to the epithelium. The volume of fiuid in the 
near-field is very greatly reduced because the cilium is much closer to the epithelium, and hence the distance from the 
image system is comparable to r even for relatively small values of r. 

In the inner and near-fields, particles follow the rotational cilium motion; the near-field however shrinks during the 
rotational beat, so a significant volume of fiuid is subject to a relatively strong 0(l/r) near-field flow during the effective 
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stroke, but a much weaker 0(l/r^) far- field flow during the return stroke. This leads to a 'loopy drift' of particles towards 
the left. Further out, particles are in the stresslet far- fields of both the effective and recovery strokes, however the former 
has a larger magnitude due to the stresslet strength being proportional to height above the surface. This leads to a 
time-averaged stresslet far- field with radial streamlines, as evident in simulation results [26l [21]. The predictions of a 
lower vortical flow and an upper directional flow were subsequently confirmed by observations of the flow in Kupffer's 
vesicle of zebrafish embryos ( [28] , figure u\ . 

An alternate perspective is to represent the time-averaged rotational motion by rotlets with axis of rotation tilted 
towards the posterior [20]. Replacing the conical rotation by a line distribution of rotlet singularities (figure [sj:), the 
inner-field has strength 0(l/r), as determined by integrating a rotlet distribution, the near- field has strength 0(l/r^), 
asymptotic to the strength of a rotlet, and the far-field, where the image singularities become equally important has 
strength 0(l/r^). Indeed, only the component of the rotlet parallel to the boundary will have a significant far-field effect, 
since a rotlet perpendicular to the boundary decays with 0(l/r^), as shown in figure [4[c, d). In the far-field, the rotlet 
has stresslet character, reconciling the two approaches for modelling the rotational cilium movement. 

Figure [g] from Hirokawa et al. [11 summarises these findings: the posterior tilt of the nodal cilium results in a reduced 
flow rate during the recovery stroke compared with the effective stroke, which generates an overall flow. 

Our analysis does not however take into account the effect of the upper membrane (figure If), which will produce 
a 'back pressure' gradient enforcing mass conservation, as first considered theoretically by Cartwright et al. [20]. In 
the next section we review the regularized stokeslet boundary integral method, and report new computational results 
applying this technique to the enclosed nodal cavity. 

5 The regularized stokeslet method 

5.1 Integral equation 

For convenience in carrying out boundary integral calculations, and to give a regular flow field throughout the domain, 
including the regions containing singularity distributions, Cortez [29 introduced the 'regularized stokeslet'. This is defined 
as the exact solution to the Stokes flow equations with smoothed point-forces, 

= -Vp + /iV^n + /V^e(a^-O. V-n = 0. (15) 

The symbol ipe{x — y) denotes a cutoff- function or 'blob' with regularisation parameter e, satisfying Jj^g ip^{x)dVx = 1- 
Cortez et al. [44] showed that with the choice i/j^ix — ^) := 15e^/(87r/irJ), the regularized Stokeslet velocity tensor is 
given by 



Here and in the rest of the paper, we use the compact notation r^ = \/t^ + e^. We follow Cortez and co-authors in using 
the solution S\a^ and its counterpart near a no-slip boundary B\- given in equation (17), as the basis for our computational 
study. 

Ainley et al. [45 derived the equivalent regularized image system, which we denote B\- for flow near a no-slip 
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boundary, which we rewrite in index notation: 



BU^,^) = ^r- 



1 fdij{r^ + 2e^)+rirj 5ij(R^ + 2e^) + RiR 



'3 



Sttii V rf Rl 

^„,. d fhRi Si3{R^ + 2e^) + RiR3 

Qhe^ 



A-KhSikMR) 



dRk \ R^ i?? 

^5 (fei?,- - %i?3) j . (17) 

The tensor Ajj^ is as defined for equation dsl; R'^ := R'^ + e^ and the term ^e(^) •= 3e^/(47ri?f) is a more slowly-decaying 
blob than i/j^^ generating the regularised source dipole image. 

Cortez et al. [44 derived the equivalent Lorentz reciprocal relation, and hence a boundary integral equation. This 
is the basis for the 'regularized stokeslet method' (RSM). Again using linearity of the Stokes fiow equations and adding 
together the boundary integral for a surface S of regularized stokeslets and an array of slender body integrals for the cilia, 
we have the following equation for the fiuid velocity at location cc, 

u,{x,t) = // 5|^.(^, 0^,(1, t)d5^+V / G,j{x,^^^\s,t))f^'^\s,t)ds. (18) 

JJs ^^1 Jo 

The unknowns <l>j(^,t) and fj\s^t) denote the j-components of the stress on the membrane S at ^ e S and force per 



unit length at arclength s on the TTith cilium respectively, these being unknown a priori. As defined in section [431 ^^^ 
kernel Gij is a combination of stokeslets and quadratically- weighted source dipoles chosen to give optimal accuracy on 
the cilia surfaces. 

5.2 Numerical implementation 

As described previously [46 , we implement a constant-element method where the numerical quadrature of the kernel and 
discretisation of the unknowns are 'decoupled'. The membrane is decomposed into elements 6'[1], . . . , 6'[A/'5'], on which the 
stress is approximated by *^j[l], • • • , ^j[Ns]; the TTith cilium is decomposed into elements /^^^[l], . . . , /^^^[A/'c] on which 
the force per unit length is approximated by /j[l], . . . , fj[Nc]- The discrete system at time t is then given by 

Ns pp M Nc 

u,{x,t) = V^,M // BUx,^)dS^^y^yjj[q] I G,j{x,e'^^{s,t))ds, (19) 






where ^q is the approximate centroid of the surface element S[q] and the arclength Sq = {q — 1/2) /Nc- This is then solved 
by collocation, applying equation (19) at a? = ic^, for / = 1, . . . , (M x Nc + Ns) locations at the midpoint of each cilium 



and the approximate centroid of each surface element, for more details, see [46]. 

The cilia are discretised with intervals I[q] = {{q — 1)/Ns^q/Ns); Ns = 12 constant-force density elements were used 
per cilium, and a membrane surface mesh with 6 x 10 x 10 = 600 constant-stress elements, with geometry based on 



the projected cube mesh used by Cortez et al. ^44^ rearranged to a hemisphere and deformed using equation (20). The 
regularisation parameter for the surface mesh was taken as e = O.OIA, where A is the notional hemisphere radius, based 
on numerical tests reported previously [44l[46]. Surface integrals were performed using 12 x 12 or 4 x 4 Gauss-Legendre 
quadrature, the more refined rule being used if the evaluation point was within a distance 1.5 of the element centre. Cilia 
stokeslet integrals were taken using standard analytic formulae (see for example [46 ); weighted source dipole integrals 
were taken using Gauss-Legendre quadrature. For the latter, quadrature with 12 points was used if the distance from 

13 



the element centre was less than 0.2, quadrature with 4 points was used if the distance was between 0.2 and 1.0, and the 
integrals were neglected at greater distances, due to their rapid O {{a/r)^) decay. 

For an array of 25 cilia beating in synchrony, as depicted in figure 8 the discrete approximations fj[q] and ^j[i/]for 
j = 1,2,3, cilium number m = 1, . . . , 25, cilium interval q — 1, . . .712 and surface element v = 1, . . . , 600 results in 
3 X (600 + 25 X 12) = 2700 scalar degrees of freedom at each timestep. Due to the absence of any explicit time-dependence 
in the Stokes flow equations, there was no coupling between timesteps; the calculations at each timestep were independent. 
The linear system was solved by L/7-decomposition, as described in [47, Chap 2.3]. Moreover, it is only necessary to 
compute the force density and stress distributions for a single beat cycle due to periodicity. Computations were performed 
using the BlueBEAR and Babbagel clusters at the University of Birmingham. 

Once the discrete approximations to j ]^ \q\ and $j \y\ have been calculated, the velocity field at any point in the fluid 



can then be calculated using equation (19). To give a regular velocity field throughout the flow domain, the cilia slender 
body integrals being taken using regularized rather than singular kernels, the regularisation parameter e being taken to 
be the minimum value of the radius a(^s) on each segment. The accuracy of this approximation was verified by computing 



the cilium surface velocity, as described in section |4.3[ Particle tracking was performed using a predictor-corrector 2nd 
order algorithm, as described previously [26^ . 

6 Results 

As in figure [2J all results are plotted in the conventional manner, with the 'left' of the embryo on the right hand side of 
the figure, in the direction of the positive x\ axis, and flow to the 'left' of the embryo will be referred to as positive. The 
X2 axis is then posterior to anterior and the x^ axis is dorsal to ventral. 

6.1 Nondimensionalisation 

Results are presented using the following scalings: length is scaled with respect to cilium length, which in the mouse 
embryo is typically 3-5 /im. Time is scaled with respect to beat period 27r/a;, typically 0.1 s for a beat frequency of 10 Hz. 
Velocity is scaled accordingly, giving a characteristic velocity scale of 30-50 /im/s. The dimensionless cilium tip speed is 
therefore 27rsin'0, typically around 130-220 /im/s. 

These scalings entail that the correct scaling for the force per unit length fj is \igL^ and for stress ^j is \ig^ where \i 
is dynamic viscosity. Viscometry data on the extraembryonic liquid do not appear to have been presented; experiments 
are typically carried out using saline media with added serum, the viscosity being within an order of magnitude of that 
of water, so /i can be estimated as 0.001 Pa s. 

6.2 Geometry 

The upper membrane surface (xi,X2,X3) = IX^Y^Z) was defined as a hemisphere, deformed using the following trans- 
formation: 

X = x(i-|), 



Y = Y, 

\z'l-\ (20) 



3 71/2 



4 
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where (X, F, Z) are the coordinates of a hemispherical surface X^ + F^ + Z^ = A^, Z > 0, with radius A = 6, as shown 
in figure [SJ The cilium beat cycle was defined by a tilted conical rotation, as in equation (12). All results were calculated 



with tilt angle = 35° and semi-cone angle tjj = 45°, to correspond with the particle tracking study in [21^ . 

6.3 Instantaneous flow with and without the overlying membrane 

Figure |9] shows instantaneous velocity field plots in the section X2 = —1.95, the nodal geometry being as shown in figure [8J 
The fiow is shown at the apices of the recovery (a) and effective (b) strokes, that is at the highest and lowest points of 
the cilium beat cycle. The fiow is in the direction of the cilia motion in the region xs < 1.75, there being a return fiow 
in the region xs > 1.75 induced by the presence of the membrane. Due to the cilia being closer to the membrane during 
the recovery stroke, the induced fiow magnitude is smaller, producing an overall positive 'leftward' fiow close to the cilia 
and a negative 'rightward' fiow in the upper region. No instantaneous counter fiow is evident close to the no-slip plane 
Xs = 0. Figures loFc, d) show the corresponding results across the same domain but with no overlying membrane. The 
principal difference is the complete absence of the upper return fiow. Additionally, when the membrane is neglected, a 
'stresslet' type fiuid fiow is evident, with 'upstream' fiuid being drawn towards the epithelium and 'downstream' fiuid 
being pushed away from the epithelium, as reported in earlier work [2JJ. 

6.4 Mean flow 

We now consider the time- averaged fiow produced over an entire beat cycle; results were computed by averaging a beat 
cycle divided into 60 discrete intervals and calculating the velocity component ui on the sagittal midplane X2 = 0. The 



right-to- left component of the velocity field is shown in figure 10 as a shaded plot, the upper panel giving an in-situ view 
of the sagittal midplane, the lower panel showing a projection, with contours depicting the zero mean fiow files separating 
leftward (positive) and rightward (negative) mean fiow. The region of leftward fiow is smaller in size than the region of 
rightward fiow; additionally this section reveals 'pockets' of mean rightward fiow close to the epithelium, induced by the 
return stroke of the cilia. The sagittal section does not however reveal the extent of these regions; for this we examine 



the constant height xs = 0.1125 section shown in figure 11 This figure shows that throughout the majority of the node, 
the negative fiow regions are localised to the cilia; between cilia there is a multiply-connected region of relatively slower 
positive fiow. In contrast, the cilia close to the posterior edge of the node create a continuous region of negative fiow; 
additionally there is a weak negative fiow towards the left, right and anterior peripheries. The effect of this mean fiow 
profile on particles within the node is unclear; we now examine particle tracking results to give insight into this. 

6.5 Particle Transport 



Figures [T2]-[T6] show the results of particle tracking simulations, a simple model for the advection of finite-sized NVPs by 
the nodal fiow. Results are shown for a period of 20000 timesteps, or 333 beat cycles, corresponding to approximately 



30 seconds at 10 Hz. Figure 12 shows tracking results with initial particle position chosen for comparison with results 
given previously [21]. The particles are initially at the locations indicated by arrows; the lowest particle is at xs = 0.1, 
the highest at xs = 1.1. All particles are initially swept to the left, even those starting in the region with negative mean 



fiow (figure 11). This is because the rotational inner-field propels the particles to the region where the leftward fiow 
dominates. Once transported to the left, all particles considered were subsequently lifted above the cilia array into the 
return fiow region; indeed the particle starting closest to the epithelium was lifted the highest. Particles were generally 
not transported beyond the extent of the cilia array, and indeed particles starting higher up travelled a shorter distance 
before being returned. The eventual location of the particle varies considerably with a small change to initial position, 
as found previously [26 and observed experimentally [28] , 
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The behaviour of particles released in the extremities of the node is shown in figures 13-16 Particles released to 



the left of the node (figure 13) are generally lifted up into the return flow, although for particles released close to the 
epithelium, this process may not be complete in the time simulated (figure [l3j a)). The particle released at X3 = 1.1 is 
returned sufficiently quickly that it re-enters the positive flow region, moving around the edges of two cilia envelopes. 
Figure [M] shows the transport of particles released just to the 'right' of the cilia array. Particles are advected to the 'left', 
although the trajectories and eventual locations are highly unpredictable, with particles initially in a column, sometimes 
being temporarily trapped by one or two cilia vortices, and particles being found on opposite sides of the node after 
20000 timesteps. Figure [15] shows trajectories with initial location at the posterior edge of the node. Particles at all 
heights from X3 = 0.1 to X3 = 1.1 are found to be transported to the right, and are pushed down a short distance towards 
the epithelium. For initial height 0.5 and above, the simulations predicted particles will orbit the cilium with base at 
(—3.75, —2.25), before being carried to the left by a succession of cilia. 

Figure [16] shows that particles released in the region anterior to the cilia array are subject to very little advection over 
a period of the duration of 333 beat cycles, at heights X3 = 0.1 and X3 = 1.1; results at other heights were similar. 

6.6 Relation to experimental observations 

The instantaneous 'mainstream' flow velocity in regions between the cilia, as shown in figure [9] typically has a peak value 
of around 0.5, corresponding to 15-25 /im/s, while the particle drift speed can be two orders of magnitude slower, with 
the leftward drift phase in figure [T2]^a) having a dimensionless speed of approximately 0.09, in dimensional units this 
corresponds to 2.7-4.5 /im/s, the rightward drift phase being slower still, the dimensionless speed being around —0.036, 
corresponding to —1.1 — 1.8 /im/s. By comparison, Okada et al. [22 reported positive flow drift speeds with a mean of 
around 3.5 /im/s in the mouse and return flow drift speeds of around —1.5 /im/s; in both cases there was considerable 
range in the observations, with for example some particles having a drift speed of up to 6 /im/s. 

7 Summary and discussion: research to date and new findings 

Previous models have established a number of important features of the generation of the nodal flow: 

• The flow is generated by whirling cilia tilted towards the posterior, approximately tracing out the surface of an 
inclined cone J2Q| fT3 l [22] . 

• Fluid dynamic models predict optimal fluid flow when the axis of the cone is tilted by approximately 35° and 
the semi-cone angle is approximately 55°, in general agreement with an experimental analogue [13 and biological 
observations [13^ .22^ _28 . 

• The viscous drag produced by the cell surface results in the relatively upright effective stroke, producing a much 
larger 'near-field' where the velocity is 0(l/r) than the recovery stroke. Moreover, the strength of the 0(l/r^) 
far- field velocity is greater during the effective stroke than the recovery stroke, being proportional to the height of 
the cilium above the surface. The combined effect is for more fluid to be transported to the left by the effective 
stroke than is returned to the right by the recovery stroke, resulting in an overall leftward flow as summarised in 
the schematic figure [6] 

• The resulting particle paths are a 'loopy drift' close to the cilia, and a 'radial drift' further away. 

• From a time-averaged perspective, image singularities convert a near-field rotational motion, characterised by an 
antisymmetric rotlet tensor, into a far-field straining motion, characterised by a symmetric stresslet tensor. As 
particles are transported by the array of cilia they may move from the far-field of one cilium to the near-field of 
another and vice versa ('multiple-cilia effect'). 
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• These effects can be understood, and modelled mathematically, through the use of image singularities located 
beneath a plane boundary, as first derived by Blake ^ . 

• The multiple-cilia effect may create chaotic advection, as observed experimentally [28]. 

• The closure of the nodal cavity will create a return flow in the upper region of the node [25] . 

This study has developed the model further by combining time-dependent cilia-driven flow with a closed nodal cavity. It 
was found that: 



• 



• 



The region of negative flow is larger than the region of positive flow, covering most of the region above the cilia 
tips, and in addition for the model considered the extreme peripheries of the node beyond the cilia array. 

The presence of the membrane significantly reduces the flow rate in the ciliated region. 

Close to the epithelium, there is generally a slow positive mean flow, with a negative mean flow in enclosed regions 
around the cilia. There may also be a layer of continuous negative flow close to the epithelium at the posterior end 
of the node. 

• Particles released between the epithelium and the height of the cilia tips are generally advected to the left edge 
of the ciliated array, before being lifted into the return flow region. The trajectories taken may however be very 
unpredictable, involving detours around cilia vortices for one or many rotations, and considerable vertical mixing 
of the flow. 

• Particles released in the posterior negative flow layer are transported to the right of the node, before entering a 
leftward trajectory. This behaviour is not present in models which neglect the upper membrane. 

• Particles released at the anterior extremity of the node were not predicted to travel a significant distance over a 
simulated period of 330 beat cycles, corresponding to about 30 seconds. In this region especially, diffusion may be 
the dominant transport mechanism. 

Our model does not yet explain how NVPs may be trapped and/or broken up by cilia at the left of the node: once 
the particles have reached the edge of the cilia array they typically are transferred to the return flow. It may be that 
cilia must extend to the edge of the node in order for this to occur, and perhaps that particle break-up and endocytosis 
[38] or some other form of morphogen inactivation [20] must prevent morphogens being returned to the right and hence 
spread equally. 

Chaotic low Reynolds number advection has been predicted to assist with the feeding of sessile micro-organisms on 
the basis of a stokeslet-image system model, using similar techniques to those described in this paper [49]. While we have 
reported qualitatively 'chaotic' behaviour, it remains to quantify mathematically the degree and extent of chaos in nodal 
flow, for example through Lyapunov exponents [49 . The apparent presence of chaotic transport is perhaps surprising in 
a structure that has evolved to break symmetry in an organised way. It is unclear to the present authors to what extent 
chaos may help, hinder, or have no significant effect on the symmetry-breaking process, and how the biological system 
may have evolved either to minimise chaos, or to exploit its potential to spread particles across the left side of the node. 

8 Future directions 

The ventral node structure and equivalent organising structures vary considerably in shape, size and other detailed aspects 
across different classes of vertebrates. Fluid mechanical modelling has a role to play in determining the importance of 
different aspects of node morphology and scale. The biological system differs from our idealised mathematical model in a 
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number of respects, not least the heterogeneity in ciha beat pattern, ciUa orientation and sometimes direction, in addition 
to irregular spacing, as shown in biological data. The epithelium itself is also not a smooth planar surface, but rather 
is made up of a 'cobbled' surface of curved cells. The regularized stokeslet method described in this study provides a 
basis for models which could take this complex geometry into account. The development and growth of the node and 
expression of cilia is also a dynamic process, and the detailed effects of these dynamic changes is yet to be uncovered. 
An additional important feature not considered in this study is the effect of synchronised cilia, versus cilia beating with 
different phase, or indeed different frequencies. Results reported previously [26] suggest that this may significantly alter 
the trajectories of particles released near the cilia envelopes, although it does not significantly change the overall particle 
trajectories. Further investigation of this effect within an enclosed nodal cavity is warranted. 

The end-point of these models is the determination of how the flow may establish a morphogen gradient, and how this 
will induce asymmetric signalling. Initial mathematical and laboratory models of this phenomenon have been developed by 
a number of groups [25l [20l |22l |28] , leading to a number of important insights. In 1970, Crick [10] proposed that diffusion 
through cells, combined with a 'source-sink' process, could produce spatial gradients of morphogens in embryogenesis. 
Recently, Yu et al. [48] demonstrated that a morphogen protein, Fgf8, can form gradients by diffusing into the extracellular 
fluid, the 'sink' function being performed by the endocytosis. These experiments were performed on living zebrafish 
embryos during gastrulation, the developmental stage at which left-right symmetry breaking is initiated. The diffusion 
measurements were achieved through the use of fluorescence correlation spectroscopy; future studies with novel techniques 
such as this will allow accurate quantification of morphogen gradients in the node associated with left-right symmetry 
breaking. Fluid mechanics plays an important role in many aspects of both animal and plant developmental biology, 
as recently reviewed by Cartwright and co-authors [50], and presents many unexplored problems at the interface of 
mathematical modelling and experimental biology. 

The past decade has seen significant advances in our biological understanding of the early stages of this process. These 
advances have inspired fluid dynamics studies, and on certain occasions have been anticipated by them. We look forward 
to future multidisciplinary work on this exciting and fundamental scientific problem. 
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Figure 5: The zones of influence of the cilium during the phases of its motion. The three fields, 'inner', 'near' and 'outer' 
decay differently with distance r from the cilium. The inner- field corresponds to r <C L, the near-field r ^ L and far- field 
by r ^ L. (a, b) from the perspective of the instantaneous cilium motion, the relatively upright effective stroke (a) 
entrains a larger region of fluid than the relatively less upright recovery stroke (b). The inner- field is essentially given by 
a line distribution of stokelets and source-dipoles (O(logr)), the near- field by the stokeslet itself (0(l/r)) and the far- field 
by the stokes-doublet image system (O (l/r^)) shown in figure H (c) From the perspective of a time-averaged rotation, 
the inner-field is given by a line integral of rotlets (0(l/r)), the near- field by the rotlet itself (O (l/r^)) and the far- field 
by the stresslet image system (O (l/r^)) shown in figure El The far-field velocity in all three cases is an O (l/r^) stresslet 
field. 

Figure 6: A depiction of the generation of the nodal flow (copyright figure not included in arxiv version). Cilia rotate 
anticlockwise, and due to their position on the curved cell surface, their rotational axis is tilted towards the posterior. 
This results in there being an effective leftward 'propulsive' motion where the cilium projects above the cell surfaces, 
and a rightward 'return stroke' where the cilium is close to the cell surfaces. The viscous drag associated with the latter 
reduces the effect of the cilium during its return stroke, creating a flow from right to left. 



Figure 7: A laser ablation particle tracking study of Kupffer's vesicle in zebrafish embryos proved the existence of the 
lower vortical flow and upper directional flow, in addition to chaotic particle transport (not shown) (copyright figure not 
included in arxiv version). 
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a 'Top' view, with mesh 



b 'Side' view, mesh outhne only 
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Figure 8: The geometry of the node used for the computations. Ellipses show the trajectory of the cilium tip over its 
rotational beat; the rod-shaped cilia are shown in their 'leftmost' position. For simulations in this paper, all cilia were 
synchronised, rotating in the clockwise direction viewed from above, and had the same geometric specification: tilt angle 
= 35°, semi-cone angle tjj = 45° and tilt direction cj) = 180°, corresponding to a completely posterior tilt, (a) View 
from ventral to dorsal, looking onto the node from 'above', with mesh, (b) View from 'left' to 'right', looking at the X2X3 
plane, (c) View from posteroir to anterior, looking at the ^1X3 plane. 
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a Apex of recovery stroke, with membrane 
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b Apex of effective stroke, with membrane 
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d Apex of effective stroke, no membrane 
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Figure 9: Instantaneous flow field across the plane X2 = —1.95 at the apices of (a, c) the recovery stroke and (b, d) 
the effective stroke. Conventionally, the positive xi direction is referred to as the 'leftward' direction, and appears as 
rightward in this figure, (a, b) show results for simulations with the overlying Reichert's membrane boundary, (c, d) 
show results without the membrane boundary. To interpret the velocity scale, the peak velocity vector in figure (b) has 
dimensionless magnitude 0.63, corresponding to a dimensional value of 19 /im/s for cilium length 3 /im and beat frequency 
10 Hz. By comparison, the peak velocity vector in figure (d) has magnitude 1.01, or 30; /im/s, showing that the upper 
membrane has a signficant effect even in the ciliated region. 
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Figure 10: Time- averaged 'leftward' flow component iZi, over the medial plane of the node xi = 0, shown in an isometric 
view, with mesh, and as a sagittal projection onto the X2X^ plane. Black indicates negative transport, white positive 
transport, with contours showing the zero mean flow lines 'Ui = 0, which separate the positive and negative flow regions. 
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Figure 11: Time- averaged 'leftward' flow component iZi, over the plane xs = 0.1125, close to the epithelium. 'Isometric' 
view, shows the membrane computational mesh; the 'top' view shows the projection onto the xiX2 plane, with contours 
showing the lines ui = separating regions of positive and negative mean flow. 
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Figure 12: Particle transport, initial particle position marked with arrow. The ellipses depict the envelope covered by 
the cilia tips viewed from above, (a) Initial position xi = —0.75, X2 = —3.25, xs = 0.1, showing xiX2 projection ('top') 
and xiXs projection ('front') . (b) Initial position as before, but with xs = 0.3. (c) Initial position as before but with 
Xs = 0.5. (d) Initial position as before but with xs = 1.1. 
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Figure 13: Particle transport, initial particle position marked with arrow, (a) Initial position xi = 4, X2 = 0, X3 = 0.1, 
showing X1X2 projection ('top') and xiXs projection ('front') . (b) Initial position as before, but with xs = 0.3. (c) Initial 
position as before but with xs = 0.5. (d) Initial position as before but with xs = 1.1. 
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Figure 14: Particle transport, initial particle position marked with arrow, (a) Initial position xi = —4, X2 = 0, xs =0.1, 
showing X1X2 projection ('top') and X1X3 projection ('front') . (b) Initial position as before, but with xs = 0.3. (c) Initial 
position as before but with xs = 0.5. (d) Initial position as before but with xs = 1.1. 
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Figure 15: Particle transport, initial particle position marked with arrow, (a) Initial position xi = 0, X2 = —5.5, xs = 0.1, 
showing X1X2 projection ('top') and xiXs projection ('front') . (b) Initial position as before, but with xs = 0.3. (c) Initial 
position as before but with xs = 0.5. (d) Initial position as before but with xs = 1.1. 
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Figure 16: Particle transport, initial particle position marked with arrow, (a) Initial position xi = 0, ^2 = 5, X3 = 0.1, 
showing X1X2 projection ('top') and X1X3 projection ('front') . (b) Initial position as before, but with xs = 1.1 
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